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A discrete Boltzmann model (DBM) is developed to investigate the hydrodynamic and thermodynamic non-equilibrium (TNE) 
effects in phase separation processes. The interparticle force drives changes and the gradient force, induced by gradients of 
macroscopic quantities, opposes them. In this paper, we investigate the interplay between them by providing detailed inspection 
of various non-equilibrium observables. Based on the TNE features, we define a TNE strength which roughly estimates the 
deviation amplitude from the thermodynamic equilibrium. The time evolution of the TNE intensity provides a convenient and 
efficient physical criterion to discriminate the stages of the spinodal decomposition and domain growth. Via the DBM simulation 
and this criterion, we quantitatively study the effects of latent heat and surface tension on phase separation. It is found that, the 
TNE strength attains its maximum at the end of the spinodal decomposition stage, and it decreases when the latent heat increases 
from zero. The surface tension effects are threefold, to prolong the duration of the spinodal decomposition stage, decrease the 
maximum TNE intensity, and accelerate the speed of the domain growth stage. 


1 Introduction 

Owing to the existence of complex interparticle interactions at 
the microscopic level and nonlinear interfaces between vari¬ 
ous phases/components at the macroscopic level, the hydro- 
dynamic non-equilibrium (HNE) and thermodynamic non¬ 
equilibrium (TNE) effects play a major role in shaping up the 
essential features of dynamic relaxation phenomena in multi¬ 
phase flow systems. The HNE and TNE show the features of 
the system in different aspects. The traditional Navier-Stokes 
model describes well weak HNE, but encounters difficulties 
in describing the TNE. To this purpose, a model based on the 
Boltzmann equation is preferable. 

In the past two decades, as a special discretization of the 
Boltzmann equation, the lattice Boltzmann method has car¬ 
ried substantially forward on the physical understanding of 
multiphase flows.- - — In recent studies,— - — the lattice Boltz¬ 
mann method was developed to probe the trans- and super¬ 
critical fluid behaviors or both the HNE and TNE simulta¬ 
neously in complex flows, which bring some new physical 
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insights into the fine structures in the system. Such an ex¬ 
tended lattice Boltzmann kinetic model or discrete Boltzmann 
model (DBM) should follow more strictly some necessary ki¬ 
netic moment relations of the equilibrium distribution function 
Thus, besides recovering the Navier-Stokes equation, it 
describes also the evolution of some non-conserved physical 
quantities, for example, the difference of the energies in var¬ 
ious degrees of freedom, the energy flux, etc, in the contin¬ 
uum limit. The TNE can be simply measured by the differ¬ 
ences between the non-conserved kinetic moments of fa and 
ffc, where fa is the discrete distribution function of the k -th 
group of particles. In other words, the TNE behaviors can be 
extracted dynamically from the DBM simulation without em¬ 
ploying additional analysis tools. 

In this paper, we present a DBM for multiphase flows with 
flexible density ratio, formulate new quantitative measures of 
TNE effects in the system, and utilize them to probe the phase 
separation process. 

2 DBM for thermal multiphase flows with flex¬ 
ible density ratio 

In 2007, Gonnella, Lamura and Sofonea (GLS) proposed a 
thermal lattice Boltzmann model for multiphase flows through 
introducing an appropriate interparticle force into the lattice 
Boltzmann equation.- In 2011, some of the present authors 
developed the GLS model by implementing the fast Fourier 
transform and its inverse to calculate the spatial derivatives. 11 
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As a result, the total energy conservation can be better held 
and the spurious velocities are refrained to a negligible scale 
in real simulations. In this work we further improve the DBM 
in two sides, insert a more practical equation of state and sup¬ 
plement a methodology to investigate the out-of-equilibrium 
features in the multiphase flow. 

The GLS-lattice Boltzmann equation reads as follows: 

+ Xki ' = - ~ \fki - fk?] + hi, (1) 

where hi takes the following form: 

hi = — [A + B • (¥*/ — u ) -\- (C+Cq) (\ ki — u) 2 ]/^f, (2) 

with 

A = -2(C + C*)7\ (3) 

B = •A.[ v ( J p vdw -pr)+ v - A - v (C v -")], (4) 

c = _J_{(p vdw _p7’)v-u + A: Vw—£(V-u) 2 

+^p 2 v-u+z:[-l(vp-vp)v-u 
—pVp • V(V -u) — Vp-Vu-Vp]}, (5) 

C?= 2 ^.pqpTVn (6) 

Here p, u, T are the local density, velocity, temperature, re¬ 
spectively. A = AVpVp — K(pV 2 p + |Vp| 2 /2)I — [pTVp • 
V(Ayr)]I is the contribution of density gradient to pressure 
tensor, I is the unit tensor, K is the surface tension coefficient. 
£ is the bulk viscosity. In the continuum limit, GLS model 
corresponds to, whereas is beyond, the thermohydrodynamic 
equations proposed by Onuki.— 

It is noteworthy that GLS model utilizes the van der Waals 
equation of state: P vdw = — |p 2 with fixed parameters. 

Due to numerical instabilities, the density ratio R between the 
liquid and vapor phases that the model can support is less than 
10. However, in practical engineering applications and nat¬ 
ural situations, R can vary greatly. For example, the density 
ratio of a liquid alloy system is close to 1, but that of water 
to steam is about 858 under the standard conditions. To im¬ 
prove the lattice Boltzmann model for simulating multiphase 
flows with large density ratio, extensive efforts have been 
made.— - — Among them, Yuan and Schaefer’s approach— is 
straightforward and effective. The core idea is that, through 
rearranging the effective mass in the pseudopotential model, 
more realistic equations of state, such as Redlich-Kwong,— 
Peng-Robinson, 20 and Carnahan-Starling— equation of state 
could be incorporated into the lattice Boltzmann model, which 


dramatically increased the density ratio, decreased the spu¬ 
rious currents, thereby minimizing the thermodynamic in¬ 
consistency. Similarly, in this work, through modifying the 
forcing term, i.e., replacing the term |p 2 V • u in eqn (5) by 
ap 2 V • u, the following Carnahan-Starling equation of state 
can be adopted 


P CS = PT 


1 + rj + r\ 2 - tj 3 

(i - v ) 3 


-ap 2 , 


(7) 


with 7] = bp /4, a and b are the attraction and repulsion pa¬ 
rameters. Subsequently, the total energy density becomes 
ej = pT — ap 2 + K\Vp\ 2 /2 + pu 2 /2. The Carnahan-Starling 
equation of state modified the repulsive term of van der Waals 
equation of state and obtained a more accurate representa¬ 
tion for hard sphere interactions. The incorporation of this 
equation into the DBM belongs to an improvement of phys¬ 
ical modeling. Chapman-Enskog analysis and the following 
numerical tests demonstrate that the revised model is thermo¬ 
dynamic consistent. Compared to models listed in refs. 10-12, 
our model is a thermal and compressible one that can be used 
to probe both the HNE and TNE effects. 


3 Two kinds of non-equilibrium effects 

To recover the thermohydrodynamic equations at the Navier- 
Stokes level, GLS model uses the following seven kinetic mo¬ 
ments, 

Mo q = M = P, (8) 

ki 

M?=E/^,=pu. (9) 

ki 

M^O = E \fki X ki ■ Vki = P{T+\ U • U), (10) 

ki Z Z 

M? = E/« v fa -VH = P(TI + im), (11) 

ki 

M 3 q = E-47 v « v fa' v « = p[7’(u a e i 3e y 5 i 3 r + e a u i3 e y 5 C[y 

ki 

+e a epu y 5 aj 5) + uuu], (12) 

M 3 q = L \fki y ki ■ Vfa-Vfa- = pu(27 + lu • u), (13) 

ki Z Z 

M 4 Q 2 = E \fki y ki • y kp/kpfki = P[(2T + —)7’I 

ki 

+uu(3P+^)], (14) 

where M^ n stands for that the m-th order tensor is contracted 
to a ft-th order one. Among the seven kinetic moment rela¬ 
tions, only for the first three ones, can be replaced by /&, 
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which means that in or out of the equilibrium, the mass, mo¬ 
mentum and energy conservations are kept. Replacing by 
fki in eqns (fTT1)-([T4l) will lead to the imbalance and the devia¬ 
tion 

A„=M„(/ fa -)-M^(47), (15) 

which can be used to measures the departure of the system 
from the local thermodynamic equilibrium. 

For an ideal gas system, the HNE and TNE effects are only 
induced by gradients of macroscopic quantities, also referred 
to gradient force. For multiphase flow system, the existence of 
interparticle force makes the situation a little more complex. 
The force term in the DBM equation works as the second driv¬ 
ing force. Especially, the right-hand side of eqn dTJ can be 
reorganized as 

Rhs = ~±L/fe-(1 + T d)fg] = (16) 

where 6 =-[A + B • (\ ki - u) + (C+ C q )(\ ki - u) 2 ], f^’ NEW = 
(1 + T 0)flf can be considered as a new equilibrium state 
shifted by the interparticle force. Thus, 

A F n =M n (xef k f)=M n (tI ki ) (17) 

are the non-equilibrium effects induced by the interparticle 
force, and what we measured from /& and 

An = Mnifiu) - IVCC4?) = A F n + (18) 

are the combined or the net non-equilibrium effects, where 

A? = M n (f ki ) - IVCC^f ,NEW ) (19) 

are the non-equilibrium effects induced by the gradient force. 
It is clear that, when the interparticle force disappears, the net 
non-equilibrium effects are only from the gradient force, i.e., 
A n = A^, corresponding to an ideal gas system. Note that, 
M n contain the information of u, so do A n which describe 
both the HNE and TNE effects. If we use the central moment 
MS(/*)=E/«(v fa -u)» which is only the representation of 
the thermo-fluctuations of molecules relative to u, then A* do 
not contain the effects of u, describing only the TNE effects. 

4 Simulation results and analysis 

In this section, we first validate the DBM model via two 
benchmarks; then investigate the HNE and TNE characteris¬ 
tics in both isothermal and thermal cases; finally, study the 
effects of surface tension on the thermal phase separation. 
Throughout our simulations, we set a = 2 and b = 0.4, then 
the critical density and temperature are p c = 1.30444 and 
T c = 1.88657. The fast Fourier transform scheme with 16-th 
order in precision* and the second order Runge-Kutta scheme 
are utilized to discretize the spatial and temporal derivatives, 
respectively. The model parameters are v\ = 1.0, V 2 = 2.0, 
V 3 = 3.0, V 4 = 4.0. 



Fig. 1 Comparisons of the coexistence densities predicted by the 
DBM model and Maxwell constructions. 



t 


Fig. 2 Variations of p, pu and ej for a phase-separating process. 



Fig. 3 Profiles of macroscopic quantities for the isothermal case at 
t = 0.01. 
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4.1 Verification and validation 

To evaluate if the model can reproduce the correct ther¬ 
modynamic equilibrium of the Carnahan-Starling system, 
we simulate the liquid-vapor coexistence curves at vari¬ 
ous temperatures with 128 x 1 lattice and periodic bound¬ 
ary conditions in both directions. The initial conditions 
are (p,T,u) = (p/, 1.70,0.0), if N x /4 < x < 3N x /4; else 
(p,7» = (p v , 1.70,0.0), where pi = 2.473 and p v = 0.458 
are the theoretical liquid and vapor densities at T = 1.70. 
Parameters are T = 10 -4 , Ax = Ay = 2.222 x 10 -3 , At 
2 x 10- 5 , K = 2.1x 10- 5 , £ = 0, q = -0.004. The initial tem¬ 
perature is 1.70 but drops by 0.01 when the equilibrium state 
has been achieved. Figure 1 shows the phase diagram recov¬ 
ered from the DBM simulations and Maxwell constructions. 
The two sets of results are in accordance with each other, even 
when T drops to 1.05, corresponding to R = pi/p v = 255. 
Clearly, it proves that the DBM is capable of handling mul¬ 
tiphase flows with large density ratio as well as ensuring ther¬ 
modynamic consistency. 

Figure 2 illustrates variations of the density Ap, the mo¬ 
mentum A(pu) and the total energy Aej for a thermal phase- 
separating process calculated from the DBM model. The ini¬ 
tial conditions are (p,T,u) = (1.5 + A, 1.0,0.0), where A is 
a density noise with an amplitude of 0.001. The remain¬ 
ing parameters are N x = N y = 128, Ax = Ay = 5 x 10 -3 , 
At = 5 x 10- 5 , T = 3 x 10" 3 , K = 5 x 10“ 5 , q = -0.002. It 
is observed that, even when the initial state is quenched much 
lower than the critical temperature, Ap and A(pu) maintain 
totally to machine accuracy. A er(t) fluctuates around its ini¬ 
tial value when t < 0.4, then keeps nearly to a constant. The 
initial fluctuation is due to the numerical discretization errors 
induced by the emergence of numerous interfaces during the 
spinodal decomposition stage. However, the maximum devia¬ 
tion of ej is 2 x 10 -7 , indicating that the DBM is adequate to 
guarantee energy conservation. 

4.2 Non-equilibrium characteristics: isothermal and 
thermal cases 

For simplicity, we first examine the non-equilibrium effects in 
one-dimensional isothermal case. We set the equilibrium den¬ 
sity profile at T = 1.74 (see Fig. 3) as the initial state PinitiaiM- 
When simulation starts, the system is suddenly quenched to 
T = 1.27 and fixed at this temperature during the whole pro¬ 
cedure. Physically, the isothermal results correspond to the 
case where the latent heat of phase transition approaches zero. 
Parameters are consistent with what we used in Fig. 1 except 
for t = 10 -3 and ^=l.lx 10 -5 . Profiles of the macroscopic 
quantities at t = 0.01 are exhibited in Fig. 3. We see the HNE 
as below. The decrease in temperature leads to the appearance 
of pressure gradients near the liquid-vapor interfaces which 
drive the vapor phase flows to the liquid side. As a result, the 


liquid (vapor) phase increases (decreases) its density, then the 
phase separation phenomenon takes place. 

Figure 4 displays the total non-equilibrium manifestations 
A 2 , A 3 , A 31 , and A 4 2 for Fig. 3, which suggests the following 
information during the procedure deviating from thermody¬ 
namic equilibrium: ( 1 ) due to the initial fields are symmetric 
about the vertical line x = N x / 2, the non-equilibrium behav¬ 
iors are also symmetric (for A 2 and A 4 2 ) or antisymmetric 
(for A 3 and A 31 ) about the same line; ( 2 ) the non-equilibrium 
effects are mainly around the liquid-vapor interfaces where 
the gradients of macroscopic quantities and interparticle force 
arise, and attain their maxima (minima) at the point of the 
maximum density difference <5p max . This can be interpreted 
as follows. In the first panel, according to the nature of 
ft?, we have = p(T + u 2 x ), then \ 2xx = Mjxx ~ x 
(T + u 2 )8p +pS(u 2 ). Quantitatively, (T + u 2 )8p is the lead- 
ing part of A 2 XX , then A 2 XX Sp approximately. Therefore, 
when Sp > 0 , A 2 XX > 0 , when Sp < 0 , A 2 XX < 0 (see the two 
troughs positioned at x = 32 and 96 for the vapor phases with 
decreasing densities). We also see that A 2 yy < A 2 XX . Numeri¬ 
cally, it is because A 2 ;y 3 ; <*TSp (due to u y = 0 ), and physically, 
the non-equilibrium driving force 4 / acts only along the x axis 
and induces stronger non-equilibrium effects. Behaviors of 
A 4; 2 can be analyzed in a similar way; (3) A^xxx shows a neg¬ 
ative peak and a positive one with different amplitudes in the 
left half part of the computational domain, so do A 3 JCXV and 
A 3 i x . Since A 2 XXX ©c S(pu x ), around the left interface Su x > 0 
while Sp is negative at first, but positive later owing to phase 
separation, so at first A 3 XXX < 0 and then A 3 XXX >0; (4) there 
are some zero-components, such as A 2 xy, A 3 xxy , A^yyy, A 3 ^ 3 ,, 
and A 4 ? 2 j 9 ;. Physically, A 2 xy accounts for the shear effects, 
A 3 xxy + A 3 _ y:v:y ; = 2 A 3 5 i-y associate with the energy flux in the y 
direction. Since the system, or more fundamentally /& (x,y), is 
symmetrical about the y axis without gradients of macroscopic 
quantities, there are neither shear effects nor energy flux along 
this direction. 

To probe the non-equilibrium effects induced by the inter¬ 
particle force and the gradient force, respectively, we present 
the x component of A F and A G in Fig. 5. Two remarkable fea¬ 
tures can be observed. Firstly, A F is opposite to and stronger 
than A g . Physically, the interparticle force is the active force 
which drives the system evolution and increases the gradients 
of macroscopic quantities, while gradient force is a passive 
one dissipated partly by viscosity, thermal diffusion and sur¬ 
face tension, etc. The former is stronger than the latter before 
attainment of the final thermodynamic equilibrium state. Sec¬ 
ondly, different from A, A F and A G achieve their maxima or 
minima at the middle of the interface, where the gradients of 
macroscopic quantities own the greatest values (except for u x ). 
This is correct since the forcing term is operating through gra¬ 
dients of macroscopic quantities. 

To examine the effects of latent heat, now we go to the ther- 
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Fig. 4 The total non-equilibrium manifestations A 2 , A 3 , A 3 \ and A 4 2 for the isothermal case as shown in Fig. 3 (right Y-axis). The density 
profile at t = 0.01 is also shown in each plot to guide the eyes (left Y-axis). 





Fig. 5 The x component of non-equilibrium manifestations induced by the interparticle force and the gradient force for the isothermal case at 
t = 0.01 as shown in Fig. 3. The density profile at the same time is also shown to guide the eyes. 
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mal case, shown in Fig. 6. Parameters and the initial state are 
consistent with what we used in the isothermal case. When 
simulation starts, the system is suddenly quenched to T = 1.0, 
but T(t) is free. So p{x) and u x (x) exhibit similar behaviors 
with Fig. 3, except for T(x). Due to the release (absorption) 
of latent heat, the temperature of the liquid (vapor) phase over 
the whole domain increases (decreases). This is the distinc¬ 
tive feature compared to isothermal case where latent heat is 
zero or exchanged with the connecting heat bath. Figure 7 
depicts the corresponding TNE quantity Aj$ where the central 
moment is employed. It is noteworthy that at t = 0.022, 
the system is farthest away from the thermodynamic equilib¬ 
rium and | AJ | has the largest value. The trace of asso¬ 
ciates with the internal kinetic energy. Compared to the first 
panel in Fig. 4, two prominent differences can be found. At 
first, in the thermal case .A^, but in the isothermal 

case, |A£J > \^ 2 yy\ (theoretically A2 = A2). Comparisons of 
A 2 in both cases demonstrate that the internal kinetic energy 
are anisotropic in different degrees of freedom and the energy 
equipartition theory is broken down under the non-equilibrium 
case. While, in the thermal case, the internal kinetic energy at 
each point is conserved, thereby A^ + A^ = 0. Secondly, 
the non-equilibrium effects in isothermal case are much pro¬ 
nounced since the gradient force and gradients of macroscopic 
quantities are much stronger than those in the thermal case due 
to the loss of latent heat. 

In Fig. 8, we exhibit the time evolutions of the averaged 
A^ = iL,,A£(U) (i = l,...,N x ;j = l,...,N y ), the av- 
eraged A^, and the averaged A^ (defined similarly) for the 
thermal case, which reveals the following scenarios. Under 
the combined actions of the interparticle force and gradient 
force, the TNE behaviors appear, increase quickly and arrive 
at their maxima. The interparticle force is stronger than the 
gradient force, hence | A* f | > |A* g |. These two kinds of driv¬ 
ing forces influence, compete and balance partly with each 
other, resulting in the overshoot, oscillation and decay in A*. 
Finally, when the system arrives at its steady state, i.e., the hy¬ 
drodynamic equilibrium state, the two forces balance totally 
with each other, and consequently the net TNE effects vanish, 
A* = A*° + A* f = 0. 

4.3 Effects of surface tension on thermal phase separa¬ 
tion 

It is well known that, when a system is instantaneously 
quenched from a disordered state into a coexistence one, the 
fluids undergo two TNE stages: W 1 the early spinodal de¬ 
composition stage and the later domain growth stage, then ap¬ 
proach the finial totally separated equilibrium state. Previous 
studies focus mainly on the domain growth law at the second 
stage.Due to the emergence of large variety of complex 
spatial patterns during phase separation, especially during the 



Fig. 6 Profiles of macroscopic quantities for the thermal case at 
t = 0.022. 



Fig. 7 TNE manifestations AJ for the thermal case at t = 0.022. It is 
noteworthy that, at this moment, the system is farthest away from 
the thermodynamic equilibrium and | AJ | has the largest value. The 
density profile at the same time is also shown to guide the eyes. 



Fig. 8 Time evolutions of A^, A^ and A^ xx for the thermal case. 
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Fig. 9 Density patterns at three representative times during thermal phase separation processes, t = 0.045 (the first row), t = 0.153 (the second 
row) and t = 4.0 (the third row). From left to right, the three columns correspond to cases with K = 10 -5 , 3 x 10 -5 and 6 x 10 -5 , respectively. 


spinodal decomposition stage, how to exactly distinguish the 
two stages is an open problem; aside from, the TNE behaviors 
during the whole process are barely concerned. Recently, with 
the help of Minkowski measures, we presented a numerical 
criterion for exactly separating the two stages and quantita¬ 
tively investigated the heat conduction, viscosity and Prandtl 
number effects on thermal phase separation.* 1 Here we fur¬ 
ther give out a physical criterion and investigate the effect of 
surface tension. 

To that aim, we run simulations on symmetric phase sep¬ 
arations with various surface tension coefficients K on peri¬ 
odic 512x512 domains. The initial conditions and parameters 
are consistent with those used in Fig. 2. Figure 9 shows the 
density patterns with various K at three representative times. 
From left to right, the three columns correspond to cases with 
K = 10 -5 , 3 x 10 -5 and 6 x 10 -5 , respectively. Figure 9 mani¬ 
fests that surface tension strongly affects the pattern morphol¬ 
ogy, the speed and the depth of phase separation procedure. 
More precisely, at t = 0.045, for the case with small K = 10 -5 , 
numerous mini domains with large density difference, sepa¬ 
rated by complicated interfaces, appear, suggesting that the 
procedure has already entered the final spinodal decomposi¬ 
tion stage. While for cases of larger K , the density variance is 
quite small, decreasing with increasing K. Nevertheless when 


t = 0.153, the averaged domain size and the phase separation 
depth for the three cases are almost the same; all cases pro¬ 
ceed to the domain growth stage. As time evolves further, we 
observe that the larger the K , the faster the phase separation, 
the bigger the averaged size, the fewer the number of domains 
and the wider the interface. Summarizing, Fig. 9 demonstrates 
that the surface tension effects prolong the spinodal decompo¬ 
sition stage but accelerate the domain growth stage. 

These results are further confirmed by the time history of 
the characteristic domain size R(t), l ^ n plotted in Fig. 10(a). 
The R(t) curves behave similarly and distinguish approxi¬ 
mately the phase separation process into two stages. At the 
first stage, it increases and arrives at a platform marked by 
the green arrow. In fact, the marked point corresponds to the 
end of the spinodal decomposition stage. Vladimirova et al. 22 
pointed out that the plateau depends on the depth of tempera¬ 
ture quench and the intensity of random noise. Here we find, 
it also depends on the surface tension. The larger the surface 
tension, the longer the duration t$D of the spinodal decompo¬ 
sition stage, and the larger the domain size for the spinodal 
decomposition stage Rsd- Our results are consistent with the 
theoretical analysis, neglecting heat conduction and viscous 
effects. 23 Essentially, phase separation is a process, through 
which the potential energy transforms into the thermal energy 
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Fig. 10 (a) Evolutions of the characteristic domain sizes R for the procedures shown in Fig. 9. (b) Evolutions of the boundary length L and the 
xx component of some TNE manifestations for the phase separation process with K = 10 -5 . (c) Evolutions of the boundary lengths L (solid 
curves) and the corresponding TNE intensities D (curves with solid symbols) for phase separation processes with various K. Here 1,3,6,18 
labeled on the L-curves indicate cases with K = 10 -5 , 3 x 10 -5 , 6 x 10' 5 ,..., 1.8 x 10 4 , respectively, (d) Duration of the spinodal 
decomposition stage t$D and the maximum TNE intensity D max as functions of K. 
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and the interfacial energy. Under the action of interparticle 
force, a liquid (vapor) embryo is continuously gaining (los¬ 
ing) molecules due to condensation (evaporation), then the in¬ 
terface emerges and part of the potential energy transforms 
into the interfacial energy that is proportional to K. There¬ 
fore, an increasing K means an increasing interfacial energy, 
an increasing t$o required for completing such an energy con¬ 
version process. On the other hand, the surface tension always 
resists the appearance of new interface to minimize the inter¬ 
facial energy. The larger the surface tension is, the stronger 
the resistance is and the longer it takes for sharp interfaces to 
form. 

Afterwards, in the domain growth stage, under the action of 
surface tension, small domains merge together to minimize the 
free energy, which naturally leads to the continuous growth in 
R(t). The slopes of R(t) curves, corresponding to the phase 
separation speeds during the domain growth stage udg , in¬ 
crease with K. Therefore, at the domain growth stage, the 
phase separation process is remarkably accelerated by the sur¬ 
face tension. Specifically, the R(t) curve with K = 6 x 10 -5 
crosses with the other two at t = 0.153, then rises quickly 
and exceeds the former two. When K varies from 10 -5 to 
3 x 10 -4 , the dependence of usd on K can be fitted by 

u DG = e + fK — C gK ) 2 + ( hK)\ (20) 

withe = 0.00764,/= 1.51 x 10 2 , g = 8.06 x 10 2 , h= 1.02 x 
10 3 , as shown in the legend of Fig. 10(a). Our results show 
qualitative agreement with theoretical analysis,— simulations 
by smoothed particle hydrodynamics ,- 5 and lattice Boltzmann 
study— for isothermal case. 

To accurately determine the tso, in Fig. 10(b) we monitor 
the time evolution of the second Minkowski measure: bound¬ 
ary length L(t) for the density threshold = 1.70 for which 
the density pattern has the largest boundary length. Mean¬ 
while, some TNE manifestations are exhibited in the same 
panel. It is interesting to note that the peak of the L(t) curve 
exactly coincides with the peaks or troughs of the TNE curves. 
Each nonzero component of A or A* describes the TNE from 
its own side. To roughly and averagely estimate the deviation 
amplitude from the thermodynamic equilibrium, we further 
define a “TNE strength" 

D = +Af + A; 2 1 +AJ 2 2 . ( 21 ) 


We may also use y A 2 + A 3 + A 2?1 + A 4 2 (or its F, or G com¬ 
ponent). In general, the DBM equation is dimensionless, so 
do A and D. D = 0 indicates that the system is in thermo¬ 
dynamic equilibrium and D > 0 out of the thermodynamic 
equilibrium. Shown in Fig. 10(c) are the time evolutions of 
L(t) (solid curves) and D(t) (curves with solid symbols, cal¬ 
culated from A* f ) for various K , where 1,3,6,..., 18 labeled 


on the L(t )-curves indicate cases with K = 10 -5 , 3 x 10 -5 , 
6 x 10 _ 5 ,...,1.8 x 10 -4 , respectively. To be seen is a per¬ 
fect coincidence between the peaks of L(t) and D(t) in pairs. 
Therefore, the time evolution of D(t) provides a convenient, 
efficient and physical way to divide reasonably the spinodal 
decomposition and the domain growth stages. The left (right) 
part of the peak corresponds to the spinodal decomposition 
(domain growth) stage. In addition, compared to the morpho¬ 
logical way, the extension of the current approach to three di¬ 
mensions is straightforward. 

Figure 10(c), again, demonstrates our conclusions: during 
the spinodal decomposition stage, the larger the surface ten¬ 
sion, the longer the time delay, the smaller the slope of the 
TNE curve, and the weaker the TNE intensity. For the case 
with larger K , the longer time delay and the subsequent rel¬ 
atively mild increase in D makes the t$D longer. For exam¬ 
ple, when K = 10 -5 , t$D = 0.045, while when K increases to 
1.8 x 10 -4 , tsD increases significantly to 0.35. When K varies 
in the range [ 10 -5 ,3 x 10 -4 ], the dependence of tsD on K can 
be fitted with the following form 

tsD = ci + bK , (22) 

with a = 0.066 and b — 1.51 x 10 3 , as shown in Fig. 10(d). 
Furthermore, due to the interface, the phase separation depth, 
as well as the gradient force and interparticle force achieve 
their peak values at the end of the spinodal decomposition 
stage, the TNE intensity is the strongest at this moment. Nev¬ 
ertheless, it is also found that the surface tension effects de¬ 
crease the maximum of the TNE strength D max approximately 
in the following way 

Anax = C + dK~ 05 , (23) 

with c = —0.073 and d = 3.30 x 10 -3 , as shown in Fig. 10(d). 
Physically, the Knudsen number is usually employed to clas¬ 
sify the level of TNE, which is defined as the ratio between the 
molecular mean-free-path A and a character length L at which 
macroscopic variations are of interest. For a phase separation 
process, we can take L to be roughly the domain size at the end 
of the spinodal decomposition stage, Rsd *- 7 Thus the mean 
Knudsen number Kn = A/2 Rsd- As displayed in Fig. 10(a), 
Rsd increases with K , thus, Kn and the TNE strength decrease 
with K oppositely. Numerically, a larger K will broaden the 
interface width, reduce the gradient force and refrain the TNE 
intensity. 

5 Conclusion 

An energy-conserving discrete Boltzmann model for multi¬ 
phase flow system with flexible density ratio is developed 
and utilized to study both the hydrodynamic non-equilibrium 
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and thermodynamic non-equilibrium effects in phase separa¬ 
tion processes. Efficient parallel implementation and abil¬ 
ity to capture the non-equilibrium effects are the two ad¬ 
vantages of the discrete Boltzmann model on computational 
and physical sides, respectively. Besides being helpful for 
better understanding the hydrodynamic non-equilibrium be¬ 
haviors in the phase separation process, the thermodynamic 
non-equilibrium effects permit to formulate a physical cri¬ 
terion to separately analyze the spinodal decomposition and 
domain growth stages. This work marks a preliminary step 
towards a deeper understanding of hydrodynamic and ther¬ 
modynamic non-equilibrium effects on phase-separation phe¬ 
nomena. Much scope is left for future investigations in the 
field. 
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